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A  3D  integrated  numerical  model  is  constructed  to  evaluate  the  thermal-fluid  behavior  and  thermal  stress 
characteristics  of  a  planar  anode-supported  solid  oxide  fuel  cell  (SOFC).  Effects  of  anode  porosity  on  per¬ 
formance,  temperature  gradient  and  thermal  stress  are  investigated.  Using  commercial  Star-CD  software 
with  the  es-sofc  module,  simulations  are  performed  to  obtain  the  current-voltage  (/-V)  characteristics 
of  a  fuel  cell  as  a  function  of  the  anode  porosity  and  the  temperature  distribution  within  the  fuel  cell 
under  various  operating  conditions.  The  temperature  field  is  then  imported  into  the  MARC  finite  element 
analysis  (FEA)  program  to  analyze  thermal  stresses  induced  within  the  cell.  The  numerical  results  are 
found  to  be  in  good  agreement  with  the  experimental  data.  It  is  shown  that  the  maximum  principal 
stress  within  the  positive  electrode-electrolyte-negative  electrode  (PEN)  increases  at  a  higher  current 
and  a  higher  temperature  gradient.  It  is  recommended  that  the  temperature  gradient  should  be  limited 
to  less  than  10.6°Cmm-1  to  maintain  the  structural  integrity  of  the  PEN. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

With  growing  international  concern  regarding  the  depletion  of 
the  world’s  natural  resources  and  the  need  for  greater  environ¬ 
mental  protection,  the  problem  of  developing  renewable  energy 
sources  has  attracted  significant  attention  in  recent  years.  Solid 
oxide  fuel  cells  (SOFCs),  in  which  energy  is  produced  via  an  elec¬ 
trochemical  conversion  process,  have  a  high  efficiency,  good  fuel 
flexibility,  low  emissions,  and  good  long-term  stability.  As  a  result, 
they  are  widely  regarded  as  an  ideal  solution  for  meeting  the 
future  energy  needs  of  both  industrial  and  private  consumers. 
However,  SOFCs  require  a  high  operating  temperature  in  order  to 
induce  an  internal  reforming  reaction  within  the  anode.  Conse¬ 
quently,  various  temperature-related  problems  exist  in  ensuring 
their  long-term  mechanical  integrity,  including  high  degradation 
rates,  high-temperature  gradients,  significant  thermo-stresses,  and 
so. 

As  a  result,  the  literature  contains  many  experimental  and 
theoretical  investigations  into  various  aspects  of  typical  SOFCs. 
For  example,  in  [1-4],  X-ray  reflection  and  ring-on-ring  methods 
were  used  to  measure  the  residual  stress  and  flexural  strength 
of  the  PEN  within  an  SOFC  at  room  or  elevated  temperatures. 


*  Corresponding  author  at:  Institute  of  Nuclear  Energy  Research,  No.  1 000,  Wun- 
hua  Rd.,  Jiaan  Village,  Longtan  Township,  Taoyuan  County  32546,  Taiwan  (R.O.C.). 
Tel.:  +886  23  4711400;  fax:  +886  23  4711409. 

E-mail  address:  kwang@iner.gov.tw  (L.-K.  Chiang). 

0378-7753 /$  -  see  front  matter  ©  2009  Elsevier  B.V.  All  rights  reserved, 
doi:  10.1 01 6/j.jpowsour.2009.1 0.011 


Many  researchers  have  utilized  experimental  techniques  to  exam¬ 
ine  the  effects  of  the  operating  parameters  on  the  behavior  and 
performance  of  typical  SOFCs  [5-9].  Other  researchers  have  used 
numerical  methods  to  examine  the  thermal-fluid  behavior  and 
thermal  stress  characteristics  of  different  SOFC  structures,  includ¬ 
ing  single  units  [10-12],  single  cells  [13-15],  and  single/multi-cell 
stacks  [16-18].  For  example,  Yakabe  et  al.  [10-12]  used  Star-CD, 
a  commercial  computational  fluid  dynamics  (CFD)  tool,  to  cal¬ 
culate  the  fuel  concentration,  temperature  and  current  density 
distribution  in  a  planar  SOFC.  The  temperature  field  within  the 
SOFC  was  then  input  to  a  finite  element  analysis  (FEA)  program 
(ABAQUS)  to  evaluate  the  thermal  stress  induced  with  the  SOFC 
structure.  The  results  showed  that  a  co-flow  pattern  was  advanta¬ 
geous  in  reducing  the  temperature  difference  between  the  various 
components  in  the  SOFC  structure,  thereby  minimizing  the  inter¬ 
nal  stresses  within  the  device.  Liu  et  al.  [14,15]  performed  CFD 
simulations  to  calculate  the  reaction  rate,  fuel  concentration,  gas 
flow  rate,  temperature  distribution,  and  current  density  in  a  sin¬ 
gle  cell  tested  under  various  operating  conditions.  The  resulting 
thermally  induced  stresses  were  then  analyzed  using  commercial 
FEA  software  (MARC).  Khaleel  et  al.  [17,18]  integrated  commer¬ 
cial  MARC  FEA  code  with  an  in-house  electrochemical  module 
to  accomplish  the  modeling  and  characterization  of  planar-type 
SOFCs.  Selimovic  et  al.  [19]  utilized  an  in-house  tool  to  evalu¬ 
ate  the  electrochemical  and  thermal  performance  of  a  bipolar 
planar  cell  and  substituted  the  temperature  field  data  into  the 
commercial  FEMLAB  program  to  calculate  the  corresponding  dis¬ 
tribution  of  the  thermal  stress.  Weil  and  Koppel  [20,21]  utilized 
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Nomenclature 

Ai  pre-exponential  factor,  Am-2 

cp  constant-pressure  specific  heat,  J(kg  I<)-1 

c  mean  constant-volume  specific  heat,  J(kg  K)-1 

d  thickness,  m 

D  mass  diffusivity,  m2  s-1 

[D]  elasticity  matrix 

E  Young’s  modulus,  GPa 

Fact  activation  energy,  J  mol-1 

F  Faraday  constant,  96,487  C  mol-1 

Fh  diffusional  energy  flux 

Fm  diffusional  flux 

h  static  enthalpy,  J  kg- 1 

H  heat  of  formation,  J  kg-1 

i  current  density,  A  m-2 

i'o  exchange  current  density,  A  m-2 

k  electrochemical  reaction  rate  constant 

[k]  stiffness  matrix 

M  molecular  weight,  g 

m°  initial  weight  fraction 

P  pressure,  Pa 

p  porosity  after  reduction 

R  universal  gas  constant,  8.3 1 4 J(mol  I<)-1 

Ri  resistance,  £2 

Rm  rate  of  mass  production  or  consumption 

T  temperature,  K 

U[  local  velocity,  m  s-1 

V  operating  voltage,  V 

V  gas  flow  rate,  ml  min-1 

V  mass  fraction 

Ym  mass  fraction  of  the  different  species  components 
with  molecular  weight 

Greek  symbols 
p  density,  kg  m-3 

r  stress  tensor 

[cr]  stress  vector 

[s]  strain  vector 

s  porosity 

8  point  deflection 

Subscripts 

a  anode 

c  cathode 

H2  hydrogen 

02  oxygen 

eff  effective 

/  failed  strength 

m  mth  composition  of  fuel  reactant 


commercial  ANSYS  software  to  analyze  the  thermal  strain  and 
stress  in  the  bonded  compliant  seal  (BCS)  used  to  seal  the  cell 
and  window  frame  components  in  planar  SOFCs.  The  results  con¬ 
firmed  the  effectiveness  of  the  BCS  design  in  decoupling  the 
effects  of  differential  thermal  expansion  between  the  two  com¬ 
ponents,  but  showed  that  some  bowing  took  place  within  the 
cell. 

This  study  constructs  a  3D  numerical  model  of  an  anode- 
supported  planar  SOFC  and  then  performs  a  series  of  simulations 
using  Star-CD  code  and  the  es-sofc  module  to  investigate  the 
current-voltage  (/-V)  characteristics  of  the  fuel  cell  and  to  deter¬ 
mine  the  temperature  field  over  a  range  of  anode  porosities  and 
operating  conditions.  The  temperature  data  are  then  imported  into 
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Fig.  1.  Flow  chart  showing  overall  simulation  procedure  for  thermo- 
electrochemical/thermal  stress  analysis  of  single  SOFC  cell. 


the  MARC  FEA  program  to  compute  the  corresponding  thermal 
stress  distribution.  Fig.  1  presents  a  flow  chart  of  the  overall  simu¬ 
lation  procedure.  The  validity  of  the  numerical  model  is  confirmed 
by  comparing  the  simulated  I-V  characteristics  of  the  SOFC  with 
the  experimental  observations.  Having  confirmed  the  validity  of 
the  model,  simulations  are  performed  to  investigate  the  maximum 
principal  stress  induced  within  the  PEN  as  a  function  of  the  oper¬ 
ating  voltage,  temperature  gradient  and  anode  porosity. 


2.  Experimental 

2.1.  Single-cell  performance  tests 

To  validate  the  performance  of  the  numerical  model,  a  series 
of  experiments  was  performed  using  the  single-cell  planar  SOFC 
test  station  constructed  at  the  Institute  of  Nuclear  Energy  Research 
(INER)  in  Taiwan  (see  Fig.  2(a)).  During  the  experiments,  hydrogen 
gas  was  supplied  to  the  anode  compartment  of  the  fuel  cell,  while 
air  was  supplied  to  the  cathode  compartment.  A  cell  with  dimen¬ 
sions  of  10  cm  x  10  cm  was  sandwiched  between  two  alumina  gas 
distributors  with  Pt  current  collecting  meshes.  As  shown  in  Fig.  2(a), 
the  cell  incorporated  Pt  and  Ni  meshes  on  the  cathode  and  anode 
sides,  respectively,  to  act  as  electrical  collectors  for  the  voltage  and 
current  probes  used  to  analyze  the  I-V  characteristics  of  the  fuel 
cell. 

The  oxidant  and  fuel  gases  were  flowed  separately  to  the  cen¬ 
ter  of  the  upper  and  lower  surfaces  of  the  PEN  via  circular  inlet 
pipes  and  were  then  distributed  over  the  PEN  surface.  The  effec¬ 
tive  area  of  the  PEN  was  around  81  cm2.  The  NiO/YSZ  anode  was 
gradually  reduced  via  an  oxidation  process  during  the  warming- 
up  process  performed  at  a  heating  rate  of  1  °Cmin-1  with  a  fuel 
gas  composition  of  10%H2/90%N2.  When  the  temperature  of  the 
PEN  reached  the  preset  point  of  775  °C,  the  H2  concentration  was 
gradually  increased  to  1 00%  hydrogen  and  humidified  at  room  tem¬ 
perature  by  bubbling  method  in  steps  of  10%  every  10  min.  Once 
the  electrochemical  response  of  the  cell  converged  to  steady-state 
conditions  (as  indicated  by  a  PEN  temperature  of  800  °C  and  a  H2 
concentration  of  100%),  the  I-V  characteristics  of  the  cell  were 
measured  using  the  current  and  voltage  probes  attached  to  the 
Pt  and  Ni  meshes  on  the  cathode  and  anode  sides  of  the  SOFC, 
respectively. 
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Fig.  2.  (a)  Schematic  illustration  of  cell  test  station  and  (b)  schematic  illustration  of 
3D  test  cell  model. 


2.2.  Flexural  strength  measurement 


The  high-temperature  mechanical  properties  of  the  various 
components  within  the  SOFC  play  an  essential  role  in  determin¬ 
ing  its  mechanical  integrity  under  normal  operating  conditions.  Of 
all  the  components  within  the  SOFC,  the  PEN  is  particularly  liable 
to  embrittlement  under  high  temperatures.  As  a  consequence, 
determining  the  flexural  strength  of  the  PEN  over  a  wide  range 
of  operating  conditions  is  essential  in  ensuring  the  reliability  and 
performance  of  the  device.  In  this  study,  the  flexural  strength  of 
various  reduced  and  unreduced  ASC  (anode-supported  cell)  spec¬ 
imens  was  measured  by  performing  ring-on-ring  biaxial  strength 
tests  in  accordance  with  the  test  procedure  laid  down  in  ASTM  C- 
1499.  According  to  ASTM  C-1499,  the  equi-biaxial  strength,  ay,  of  a 
specimen  is  given  by 


(Tf  = 


3  Fl 
2jrh2 


(1  -  v) 


2D2 


+  (1  +  v)  In 


Ds 

Dl 


(1) 


where  h  is  the  thickness  of  the  specimen,  FL  is  the  applied  force,  and 
D,  Dj  and  Dl  are  the  diameters  of  the  specimen,  support  ring  and 


Table  1 

Input  parameters  for  SOFC  simulations. 


PEN  thickness 

600  (Jim 

Anode  thickness 

560  (Jim 

Electrolyte  thickness 

5  fxm 

Cathode  thickness 

35  fxm 

Rib  thickness 

1  mm 

Gas  distributors  thickness 

30  mm 

Operating  temperature 

800  °C 

Operating  pressure 

1  atm 

Cell  temperature 

800  °C 

Anode  fuel 

h2 

Cathode  fuel 

Air 

Fuel  inlet  temperature  (anode  and  cathode) 

760  °C 

Pre-exponential  factor,  A, 

128 

Activation  energy,  Eact 

550J  mol 

load  ring,  respectively.  In  general,  the  Young’s  modulus  of  a  circular 
specimen  can  be  obtained  from  ASTM  C-1499 


3F(l-v2)D2  /D|  I-  (1  -  v)(D|  -  Df) 
Sn80h3  yD2  +  2(1 +y)D2 


where  <50  is  the  deflection  of  the  center  point  of  the  specimen. 
However,  in  the  ring-on-ring  tests  performed  in  this  study,  the 
deflection  is  measured  at  the  point  where  the  load  ring  contacts 
the  specimen  rather  than  at  the  center  point.  As  a  result,  Eq.  (2)  is 
modified  to  the  form  [27] 


3F(1  -  v2)D2 
SttSiH3 

-(,+21"!))’  (31 


/D|+  ( 1  -  v)(D|  -  d\ )(d|  -  d\ ) 
\D2  [  2(1  +  v)D2D2 


where  8L  is  the  deflection  at  the  point  of  contact  between  the  spec¬ 
imen  and  the  load  ring. 


3.  Numerical  method 


3  A.  Geometric  model 


In  this  study,  the  thermal-fluid  behavior  and  electrochemical 
reactions  within  the  SOFC  cell  were  analyzed  numerically  using 
a  3D  model.  The  temperature  field  computed  in  the  simulations 
was  then  imported  to  the  MARC  FEA  program  in  order  to  cal¬ 
culate  the  thermal  stress  within  the  PEN.  Fig.  2(b)  illustrates  the 
basic  components  within  the  numerical  model,  namely  a  PEN  elec¬ 
trode  and  two  alumina  gas  distributors  (one  on  either  side  of  the 
PEN).  The  gas  flow  configuration  consists  of  multiple  ribs  measuring 
3  mm  x  3  mm  x  1  mm  on  the  alumina  gas  distributors  and  a  circular 
gas  inlet  channel  in  the  center  of  both  gas  distributors.  The  major 
parameters  of  the  numerical  model  are  summarized  in  Table  1. 


3.2.  Thermal-fluid  model 


In  the  present  simulations,  the  Navier-Stokes  equation  is  solved 
for  each  of  the  computational  mesh  elements  within  the  fluid 
region.  To  simplify  the  calculations,  the  following  assumptions  are 
made:  (1)  the  fuel  and  oxidant  gases  are  ideal  gases;  (2)  the  gas 
flow  is  steady  and  laminar  everywhere  in  the  cell  channels;  and  (3) 
the  anode  and  cathode  materials  are  isotropic  porous  media.  Given 
these  assumptions,  the  governing  equations  for  the  conservation  of 
mass,  momentum,  enthalpy  and  species  are  formulated  as  follows 
[14]: 

Mass  equation 


depu[  = 

dXj 


(4) 


1898 


L.-K.  Chiang  et  al.  /  Journal  of  Power  Sources  195  (2010)  1895-1904 


where  s  is  the  porosity  of  porous  media  (equals  to  zero  in  the  flow 
channel  region)  and  ut  is  the  local  velocity  of  fuel  defined  by  Darcy 
law  [28]: 

Uj  =  — —VP  (5a) 

fi/U. 

Momentum  equation 

d  .  ,  3P  £2fjLUi 

-{epujui-Tij)  =  -eWi+—,  (5b) 

where  ty  is  the  stress  tensor  relative  to  the  velocity  gradient  and 
£2fjiUilk  is  equal  to  zero  in  the  flow  channel  region.  In  addition,  the 
density  of  the  fuel  mixture  is  given  by  the  ideal  gas  law  as 


where  Ym  is  the  mass  fraction  of  the  different  species  components 
with  molecular  weight  Mm. 

Enthalpy  equation 


9  ,  ,  c  .  9P 
( phuj  +  FhJ)  =  —  +  u, 


dx. 


dt 


9P 

dx. 


+  rii 


diij 


■  +  sh- 


(7) 


In  the  above  equations,  source  terms  sh  represent  the  heat  gen¬ 
eration  from  chemical  reaction  and  the  specific  enthalpy  h  of  the 
fluid  is  defined  as  the  sum  of  the  thermal  and  chemical  components, 
i.e. 


h  =  cpT  -  c°T°  +  UYmHm,  (8) 

where  Hm  is  the  heat  of  formation  of  constituent  m  and  cp  is  the 
mean  constant-pressure  specific  heat  at  temperature  T  with  c° 
being  the  reference  specific  heat  at  a  reference  temperature  T°. 
For  solids  and  constant-density  fluids,  the  transport  equation  for 
a  specific  internal  energy  e  is  given  by 

e  =  cT  -  c°T°  +  SYmHm,  (9) 

where  c  is  the  constant-volume  specific  heat. 

Species  equation 

9 

-^-(pUjYm  +Fmj)  =  Rm,  (10) 

where  Fm  is  the  diffusion  flux  component  and  Rm  is  the  rate  of  mass 
production  or  consumption  during  the  chemical  reaction. 


3.3.  Electrochemical  model 


The  numerical  solution  procedure  includes  various  electro¬ 
chemical  reaction  and  net  enthalpy  generation  calculations, 
including  the  combustion  of  oxygen  with  hydrogen  and  the 
gas-water-shift  reaction.  In  the  computation  process,  the  amount 
of  oxygen  at  the  cathode  side  which  subsequently  reacts  with  the 
hydrogen  at  the  anode  side  in  accordance  with  H2  + (l/2)02  ->  H20 
is  determined  by  the  specified  current  density. 

The  limiting  current  densities  for  the  chemical  reactions  in  the 
anode  and  cathode,  respectively,  are  specified  as  [6] 


2FDeff(T)PH2 

RTda 


(11) 


4  FPDeff(T)  (  P system  \ 

RTdc  y  Psystem  ~  Pq2  J 


(12) 


where  D^T)  is  the  effective  diffusion  coefficient  and  is  given  by 
Deff(T)  =  D(T)  ■  s1'5,  (13) 

while  Ph2  and  Po2  are  the  partial  pressures  of  hydrogen  and  oxy¬ 
gen,  respectively,  F  is  the  Faraday  constant,  and  da  and  dc  are  the 
thicknesses  of  the  anode  and  cathode  electrodes,  respectively. 


The  electrical  performance  of  the  SOFC  is  evaluated  from  the 
fuel  and  air  compositions.  The  cell  voltage  is  calculated  as 


V(i)  —  Vjvernst  ^ Ohmic  ^Butler -Volmer  +  V cathode  +  ^ anodes 


(14) 


where  VNemst  is  the  Nernst  or  open  circuit  potential,  Vohmic  is  the 
ohmic  potential  drop,  and  VButier-voimer »  ^ anode  and  Y cathode  the 
Butler-Volmer,  anode  and  cathode  polarizations,  respectively.  The 
various  voltage  terms  in  Eq.  (14)  are  defined  as  follows: 


i/  RT  i  (  Po2, anode  \ 

VNemst  =  4p  ln  p 

\1  o2,  cathode  J 

(15) 

Y Ohmic  =  iRh 

(16) 

VButler-Volmer  =  ^  Sinh’1  (^)  , 

(17) 

Y cathode  —  > 

(18) 

V—.gln  (,-i). 

(19) 

Note  that  in  Eq.  (17),  the  exchange  current  density  io  corre¬ 
sponds  to  the  dynamic  electron  transfer  rate  under  equilibrium 
conditions  and  is  expressed  in  accordance  with  [17]  as 

i0  =  A  -  exp  (20) 

where  the  pre-factor  and  activation  energy  Eact  are  properties 
specific  to  the  electrode-electrolyte  interface  in  the  electrochemi¬ 
cal  model.  In  the  present  study,  the  pre-factor  A/  is  estimated  from 
the  experimental  I-V  curves. 

Boundary  conditions  for  the  dependent  variables  of  the  trans¬ 
port  equations  at  the  interfaces  between  different  layers  of  the 
same  domain  are  not  required.  However,  fully  developed  flow  con¬ 
ditions  are  assumed  at  the  fuel  and  air  flow  outlets,  i.e. 


dUj  dYm 
dz  dz 


(21) 


Meanwhile,  the  boundary  conditions  at  the  interfaces  between  the 
flow  field  and  the  gas  distributors  are  specified  as 


Oj  —  Ym  —  0. 


(22) 


Similarly,  the  boundary  conditions  at  the  interfaces  between  the 
flow  field  and  the  PEN  are  given  as 


Z=ZMEA 


dUi 

dz 


z—zflowfield 


& effU-i  I  —  till 


(23a) 


Fig.  3.  Comparison  of  experimental  and  numerical  results  for  variation  of  cell  volt¬ 
age  and  power  density  with  current  density. 
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Fig.  4.  Variation  of  cell  voltage  and  power  density  with  current  density  as  function 
of  anode  porosity. 


eff  9 z 


Z=ZMEA 


dYm 

dz 


-zflowfield 


Ym  I  =  Ym  | 


(23b) 


Finally,  the  boundary  conditions  at  the  interface  between  the  elec¬ 
trolyte  and  the  electrodes  are  specified  as 


3  U[ 
dz 


=  0, 


z~z electrolyte 


3  Ym 
dz 


=  0. 


z~z  electrolyte 


(24) 


Since  the  test  cell  has  a  symmetric  geometry,  the  3D  FEA  model 
constructed  in  this  study  considered  just  one  quarter  of  the  test 
cell  configuration.  The  model  was  constructed  using  PATRAN  soft¬ 
ware  (MSC  Inc.),  and  comprised  a  total  of  17,974  sold  elements 
(Type  84)  and  24,218  nodes.  Having  constructed  the  model,  the 
CFD  computational  domain  was  meshed  using  GRIDGEN.  Each  ele¬ 
ment  within  the  mesh  was  then  taken  as  a  control  volume  in  the 
Star-CD  finite-volume  solution  procedure.  In  other  words,  the  gov¬ 
erning  equations  were  numerically  integrated  over  each  control 
volume,  where  all  of  the  related  variables  and  material  properties 
were  stored  at  the  cell  center. 


3.4.  Thermal  stress  model 


transferring  the  temperature  profiles  to  MARC,  the  temperature 
data  were  read  into  the  nodes  and  were  then  interpolated  to  the 
corresponding  points  within  each  element.  Since  the  format  of  the 
temperature  data  generated  in  the  Star-CD  simulations  was  not 
directly  compatible  with  the  data  requirements  of  MARC,  a  propri¬ 
etary  subroutine  was  programmed  in  FORTRAN  77  to  resolve  the 
incompatibility  issue.  In  performing  the  MARC  analyses,  a  “touch¬ 
bonding  constraint”  condition  was  specified  between  the  PEN  and 
its  adjacent  components.  In  addition,  a  frictional  coefficient  of  0.16 
was  specified  for  the  contact  conditions  between  all  neighboring 
components  in  the  cell  [22].  Moreover,  an  assumption  was  made 
that  the  SOFC  cell  was  simply  supported  on  its  lower  surface.  As  a 
result,  the  entire  surface  of  the  lower  interconnect  of  the  cell  was 
constrained  in  the  direction  normal  to  the  lower  plane.  In  addition, 
a  uniform  load  was  imposed  on  the  surface  of  the  upper  intercon¬ 
nect  such  that  a  close  contact  was  maintained  between  the  PEN 
and  the  upper  and  lower  (cathode  and  anode)  gas  distributors. 
Finally,  an  additional  constraint  was  applied  to  the  edge  points  on 
the  symmetric  plane  to  prevent  rigid  body  motion  of  the  entire  test 
cell. 

In  a  previous  study  by  the  current  group  [15],  the  temperature- 
dependent  elastic  properties  of  a  NiO-YSZ  anode  [23]  were 
evaluated  without  considering  the  anode  porosity.  However,  in 
practice,  the  mechanical  properties  of  Ni-based  anode  materials 
depend  strongly  on  their  porosity  [24].  Thus,  in  this  study,  the 
elastic  modulus,  Poisson’s  ratio  and  biaxial  strength  of  the  Ni-YSZ 
anode  (in  both  an  unreduced  and  a  reduced  condition)  were  mod¬ 
eled  as  a  function  of  the  anode  porosity.  Radovic  [24]  and  Selcuk 
[25,26]  presented  the  effective  mechanical  properties  of  various 
materials  of  interest  for  SOFC  applications  (e.g.  YSZ,  GCO,  TZP,  NiO- 
YSZ)  based  upon  a  theoretical  treatment  of  the  porosity  effect  (i.e. 
the  composite  sphere  method,  CSM)  and  an  empirical  treatment 
(i.e.  exponential,  nonlinear,  and  linear  equations),  respectively.  For 
NiO-YSZ,  the  empirical  equations  were  shown  to  provide  better 
results  than  the  CSM  equations  [25].  Thus,  in  the  present  simula¬ 
tions,  the  Young’s  modulus,  Poisson’s  ratio  and  biaxial  strength  of 
the  NiO-YSZ  and  Ni-YSZ  anodes  are  calculated  at  different  porosi¬ 
ties  on  the  basis  of  empirical  equations.  In  evaluating  the  effect 
of  the  reduction  process  on  the  porosity  of  the  anode,  the  porosity 
after  reduction^,  was  derived  directly  from  the  initial  (unreduced) 
porosity,  po,  in  accordance  with  [24] 


The  temperature  profiles  generated  by  the  Star-CD  and  es-sofc 
simulations  were  imported  into  the  MARC  FEA  program  to  calcu¬ 
late  the  thermal  stress  distributions  within  the  SOFC  test  cell.  In 
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Fig.  5.  Effects  of  anode  porosity  on  current  density  distribution  in  PEN  at  operating  voltage  of  0.7  V. 
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Table  2 

Biaxial  flexural  strength  test  results  for  unreduced  and  reduced  ASC-type  PENs. 


Specimen  no. 

Flexural  strength  (MPa) 

Remark 

A-l 

151.7 

Unreduced 

A-2 

135.1 

Unreduced 

A-3 

102.7 

Reduced 

A-4 

117.1 

Reduced 

where  pNi  and  pNi0  are  the  densities  of  Ni  and  NiO,  respectively,  and 
m°i0  is  the  initial  weight  fraction  of  NiO  in  the  NiO-YSZ  composite 
material. 

According  to  the  theory  of  elasticity,  the  relationship  between 
the  nodal  displacement  and  the  element  strain  can  be  expressed  as 
[15] 

{£}  =  [B]  {d}  ,  (26) 

where  [B]  is  a  matrix  related  to  the  nodal  displacement  and  element 
strain. 

The  corresponding  stress  can  be  computed  as  follows: 

{or}  =  [D]{s}  =  [ D][B]{d },  (27) 


Fig.  7.  EDS  compositional  analysis  of  fracture  surface  of  specimen  A. 

4.  Numerical  and  experimental  results 


where  [D]  is  an  elasticity  matrix  expressed  in  terms  of  the  Young’s 
modulus  and  the  Poisson’s  ratio.For  a  given  temperature  gradient, 
the  thermal  strain  vector  is  given  by 

{sth}  =  AT[  ax  oty  oiz  0  0  0]T,  (28a) 

{£}  =  {£e}  +  {sP}  +  {8th},  (28b) 

where  {£},  {£e}  and  {£p}  are  the  total,  elastic  and  plastic  strain 
vectors,  respectively. 

According  to  the  principle  of  virtual  work,  a  virtual  change  of  the 
internal  strain  energy  (internal  work,  Wa)  is  offset  by  an  identical 
change  in  the  external  work  (WF)  attributed  to  the  applied  load. 
Thus,  the  force  vector  and  stiffness  matrix  can  be  derived  as  follows: 


(F)  =  J  [B]T [D]{s}dV, 

(29) 

[X]  =  f  [B]T[D][B]dV. 

(30) 

4.1.  Thermo -fluid  analysis 

To  confirm  the  validity  of  the  3D  model  and  the  numerical  solu¬ 
tion  procedure,  the  simulation  results  obtained  for  the  variations 
of  the  cell  voltage  and  power  density  with  the  current  density  were 
compared  with  the  experimental  results  obtained  from  the  test  cell 
shown  in  Fig.  3.  The  corresponding  results  are  presented  in  Fig.  3 
for  an  operating  temperature  of  800  °C,  a  hydrogen  flow  rate  of 
335  ml  min-1  and  an  oxygen  flow  rate  of  670  ml  min-1 .  While  the 
operating  temperature  is  set  at  800  °C,  the  real  operating  temper¬ 
ature  in  the  center  of  PEN  measured  by  thermal-couple  reaches 
to  820  °C  at  Vceu  =  0.7.  From  inspection,  the  maximum  deviation  of 
cell  performance  between  the  simulation  and  experimental  data 
is  found  to  be  less  than  6%  and  the  temperature  difference  in  the 
center  of  PEN  is  5  °C  at  Vceu  =  0.7.  Thus,  the  validity  of  the  numerical 
model  is  confirmed. 

Fig.  4  illustrates  the  numerical  results  obtained  for  the  variations 
of  the  cell  voltage  and  power  density  with  the  current  density  as  a 
function  of  the  anode  porosity.  It  can  be  seen  that  the  power  den¬ 
sity  is  insensitive  to  the  anode  porosity  at  low  values  of  the  current 
density  (i.e.  <600mAcm-2).  However,  at  higher  current  densities, 
the  power  density  increases  with  an  increasing  porosity.  Specifi- 


Fig.  6.  SEM  micrograph  of  fracture  surface  of  specimen  A. 


Fig.  8.  Variation  of  maximum  principal  stress  in  PEN  with  cell  voltage  for  anode 
porosities  of  0.2,  0.3  and  0.5,  respectively. 
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Fig.  9.  Temperature  gradient  distributions  (left)  and  thermal  stress  distributions  (right)  in  PEN  at  cell  voltages  of:  (a)  0.7  V;  (b)  0.8  V;  (c)  0.9  V  (note  that  anode  porosity  =  0.5 
in  every  case). 


cally,  the  output  power  increases  by  around  10%  as  the  porosity  is 
increased  from  0.2  to  0.3,  and  by  a  further  10%  as  the  porosity  is 
increased  from  0.3  to  0.5.  This  result  is  to  be  expected  since  the  fuel 
diffusion  efficiency  within  the  PEN  improves  as  the  anode  porosity 
is  increased.  However,  the  performance  improvement  reduces  to 
less  than  1 0%  as  the  porosity  is  increased  beyond  0.5  due  to  the  over¬ 
potential  produced  by  the  corresponding  increase  in  the  electronic 
resistance. 


Fig.  5  presents  the  local  current  density  distributions  in  the  PEN 
for  anode  porosities  of  6  =  0.2,  0.3  and  0.5,  respectively,  and  a  con¬ 
stant  cell  voltage  of  Vcen  =  0.7  V.  It  can  be  seen  that  in  every  case,  the 
current  density  reduces  radially  from  the  center  of  the  PEN  toward 
the  outlet  as  a  result  of  the  reduction  in  the  concentration  caused 
by  the  electrochemical  reaction.  Moreover,  it  is  shown  that  the  cur¬ 
rent  density  in  the  channel  area  is  slightly  higher  than  that  in  the  rib 
area  since  the  fuel  diffuses  directly  into  the  electrode  in  the  channel 
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Fig.  10.  Temperature  gradient  distributions  (left)  and  thermal  stress  distributions  (right)  in  (a)  gas  distributors;  (b)  PEN  (note  that  cell  voltage  =  0.7  V  and  anode  porosity  =  0.3). 


area,  but  travels  along  a  longer  diffusion  path  to  reach  the  electrode 
in  the  rib  area. 

4.2.  Flexural  strength  measurement 

Ring-on-ring  biaxial  strength  tests  were  performed  to  measure 
the  flexural  strength  of  four  reduced  and  unreduced  ASC  PEN  spec¬ 
imens.  Table  2  presents  the  maximum  flexural  stress  in  each  of  the 


Fig.  11.  Schematic  definition  of  temperature  gradient  on  PEN  surface  (one  quarter 
model). 


four  specimens.  It  is  observed  that  the  flexural  stress  in  specimen 
A  (unreduced)  decreases  by  around  25-30%  in  the  reduced  spec¬ 
imens.  Radovic  and  Lara-Curzio  [24]  showed  that  the  porosity  of 
NiO-YSZ  increases  from  20  to  38.2%  following  the  reduction  pro¬ 
cess.  In  addition,  it  was  shown  that  the  change  in  porosity  reduced 
the  maximum  flexural  strength  of  the  electrode  from  97  MPa  in  the 
unreduced  state  to  68  MPa  in  the  reduced  state.  In  other  words, 
the  reduction  process  reduced  the  fracture  strength  of  the  Ni-YSZ 
anode  by  around  30%.  This  finding  is  in  good  agreement  with  the 
experimental  results  presented  in  Table  2. 

The  fracture  surfaces  of  the  four  ASC  PENs  were  examined  by 
scanning  electron  microscopy  (SEM)  and  energy  dispersive  spec¬ 
troscopy  (EDS).  Fig.  6  shows  the  fracture  surface  of  specimen  A  (i.e. 
an  unreduced  PEN).  It  can  be  seen  that  the  fracture  surface  contains 
a  relatively  uniform  distribution  of  small  cavities.  These  cavities 
are  essential  for  gas  diffusion,  but  should  be  properly  controlled  in 
order  to  maintain  the  integrity  of  the  PEM  structure.  Fig.  7  presents 
the  EDS  analysis  results  for  the  fracture  surface  of  specimen  A,  and 
shows  that  the  dominant  materials  of  the  anode,  electrolyte,  and 
cathode  are  Ni,  Zr  and  LSM,  respectively. 

4.3.  Thermal  stress  analysis 

Tables  3a-3c  show  the  variations  of  the  maximum  principal 
stress  and  maximum  temperature  with  the  cell  voltage  (0.7-0.9  V) 
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Table  3a 

Maximum  principal  stresses  and  temperatures  at  different  cell  voltages  (poros¬ 
ity  =0.2). 


Voltage  (V) 

PEN  (MPa) 

Gas  distributors 
(MPa) 

PEN 

Tnax  (°C) 

PEN 

Armax  (°C) 

0.7 

20.7 

42.9 

821.2 

12.9 

0.8 

19.6 

40.5 

809.9 

12.8 

0.9 

18.1 

37.3 

808.3 

12.5 

Table  3b 

Maximum  principal  stresses  and  temperatures  at  different  cell  voltages  (poros¬ 
ity  =0.3). 


Voltage  (V) 

PEN  (MPa) 

Gas  distributors 
(MPa) 

PEN 

Tnax  (°C) 

PEN 

ATmax  (°C) 

0.7 

18.7 

42.9 

823.0 

12.8 

0.8 

18.4 

42.1 

811.3 

12.6 

0.9 

16.1 

36.8 

806.8 

11.9 

Table  3c 

Maximum  principal  stresses  and  temperatures  at  different  cell  voltages  (poros¬ 
ity  =0.5). 


Voltage  (V) 

PEN  (MPa) 

Gas  distributors 
(MPa) 

PEN 

Tnax  (°C) 

PEN 

ATmax  (°C) 

0.7 

15.1 

42.7 

824.9 

12 

0.8 

15.3 

43.5 

812.7 

13 

0.9 

13.1 

36.2 

805.9 

11.5 

for  anode  porosities  of  £  =  0.2,  0.3  and  0.5,  respectively.  Note  that 
the  gas  inlet  temperature  is  Tj  =  760oC  in  every  case.  In  theory, 
a  lower  operating  voltage  yields  a  higher  current  density  and  a 
higher  power  density,  and  therefore  results  in  a  higher  tempera¬ 
ture  field  and  a  greater  temperature  gradient  (AT).  Accordingly, 
Tables  3a  and  3b  show  that  the  thermal  stress  increases  with  a 
reducing  operating  voltage  for  anode  porosities  of  £  =  0.2  and  s  =  0.3, 
respectively  (see  also  Fig.  8).  (Note  that  these  results  are  consistent 
with  those  presented  by  Selimovic  et  al.  [19].)  However,  as  shown 
in  Table  3c,  for  a  higher  porosity  of  £  =  0.5,  the  maximum  thermal 
stress  is  obtained  at  the  intermediate  voltage  value  of  0.8  V  (see  also 
Fig.  8).  This  result  is  clarified  in  Fig.  9,  which  shows  that  for  a  poros¬ 
ity  of  £  =  0.5,  the  temperature  gradient  is  enhanced  at  a  voltage  of 
0.8  V  and  therefore  induces  a  greater  thermal  stress  within  the  PEN. 
In  [24],  it  was  shown  that  the  reduction  in  the  elastic  modulus  of 
the  NiO-YSZ  anode  relative  to  that  of  the  Ni-YSZ  anode  was  directly 
attributable  to  an  increased  anode  porosity.  Hence,  the  reduction 
in  the  principal  stress  with  an  increasing  anode  porosity  shown  in 


Table  4 

Maximum  principal  stresses  in  PEN  and  gas  distributors  at  different  temperature 
gradients. 


Max.  temperature  gradients  (°C)  (Ti-T4) 

PEN  (MPa) 

Gas  distributors 
(MPa) 

50 

60.1 

90.7 

75 

91.6 

137.3 

100 

123 

184.0 

Table  5 

Maximum  principal  stresses  in  PEN  and  interconnects  at  different  loads  (V=0.7, 
porosity  =  0.3). 


Load  (kg cm-2) 

PEN  (MPa) 

Interconnect  (MPa) 

0.5 

18.7  (-10.5) 

42.9  (-4.7) 

1 

18.7  (-11.1) 

42.9  (-5.0) 

5 

18.1  (-15.2) 

42.3  (-7.5) 

Fig.  8  can  be  attributed  to  lower  values  of  the  elastic  modulus,  shear 
modulus  and  Poisson’s  ratio,  respectively.  As  shown  in  Fig.  9(a)-(c) 
(right),  the  maximum  principal  stress  within  the  PEN  structure 
occurs  at  the  rib  edges  due  to  CTE  mismatches  between  the  PEN 
and  the  interconnect  and  the  effects  of  the  contact  force.  Moreover, 
the  maximum  principal  stress  in  the  gas  distributors  is  induced  in 
the  air  inlet  region  due  to  the  higher  temperature  gradient  (see 
Fig.  10(a)). 

Fig.  1 1  presents  a  schematic  illustration  of  the  temperature  gra¬ 
dient  induced  within  the  PEN.  Table  4  indicates  the  maximum 
principal  stresses  in  the  PEN  and  gas  distributors  for  tempera¬ 
ture  gradients  {T^-T4}  see  Fig.  11)  of  50,  75  and  100°C  (7.07,  10.6 
and  14.1°Cmm_1),  respectively.  Note  that  T4  is  assumed  to  be 
700  °C  in  every  case.  The  results  clearly  show  that  the  maximum 
principal  stress  increases  with  an  increasing  temperature  gradi¬ 
ent.  According  to  the  biaxial  test  results  presented  in  Table  2,  the 
average  strength  of  an  unreduced  anode-supported  PEN  is  approx¬ 
imately  135-151  MPa,  depending  on  its  composition.  However, 
the  average  strength  reduces  to  approximately  102-1 17  MPa  after 
the  reduction  process.  Table  4  shows  that  the  maximum  princi¬ 
pal  stress  induced  in  the  PEN  by  a  temperature  gradient  of  75  °C 
(10.6°Cmm_1)  is  91.6  MPa.  This  value  is  close  to  the  minimum 
fracture  strength  of  the  PEN  (i.e.  102  MPa,  see  Table  2).  Thus,  in 
designing  and  operating  the  SOFC  cell,  the  temperature  gradient 
on  the  PEN  surface  should  be  carefully  controlled  such  that  it  does 
not  exceed  75  °C  (10.6  °C  mm-1 )  in  order  to  ensure  the  mechanical 
integrity  of  the  device. 

Table  5indicates  the  maximum  principal  stress  within  the 
PEN  and  the  gas  distributors  for  uniform  loads  of  0.5,  1.0,  and 
5.0  kg  cm-2,  respectively.  As  expected,  the  results  show  that  the 
application  of  a  uniform  load  to  the  upper  interconnect  induces  a 
compressive  stress  within  both  the  PEN  and  the  gas  distributors. 
However,  it  can  be  seen  that  the  magnitude  of  the  compressive 
stress  increases  only  slightly  with  an  increasing  load.  Thus,  it  is 
inferred  that  the  stress  induced  within  the  SOFC  cell  is  the  result 
primarily  of  the  temperature  gradient  (Ti  -T4)  and  the  effects  of  CTE 
mismatches  between  the  various  components  of  the  cell. 

5.  Conclusion 

This  study  has  performed  a  series  of  numerical  simulations  to 
characterize  the  thermal-fluid  behavior,  electrochemical  reactions 
and  thermal  stress  properties  of  a  planar  anode-supported  SOFC  for 
various  values  of  the  operating  voltage,  temperature  gradient  and 
anode  porosity.  The  validity  of  the  numerical  model  has  been  con¬ 
firmed  by  comparing  the  simulated  I-V  characteristics  of  the  fuel 
cell  with  the  experimental  results.  The  numerical  results  support 
the  following  major  conclusions: 

1 .  For  anode  porosities  in  the  range  0.2  <  e  <  0.5,  the  cell  voltage  and 
power  density  increase  with  increasing  porosity. 

2.  A  lower  operating  voltage  increases  the  local  current  density  and 
temperature  gradient  and  therefore  induces  a  higher  principal 
stress  within  the  PEN. 

3.  The  application  of  a  uniform  load  to  the  upper  interconnect  in 
the  SOFC  device  induces  a  compressive  stress  within  the  PEN. 
The  maximum  principal  stress  increases  only  slightly  with  an 
increasing  load.  Thus,  it  is  inferred  that  the  stress  within  the 
PEN  is  dominated  by  the  effects  of  the  temperature  gradient  and 
CTE  mismatches  between  the  various  components  of  the  cell, 
respectively. 

4.  The  thermal  stress  within  the  PEN  is  smaller  at  a  lower  tem¬ 
perature  gradient.  In  practice,  temperature  gradient  of  the  PEN 
should  be  less  than  1 0.6  °C  mm-1  in  order  to  guarantee  the  struc¬ 
tural  integrity  of  the  SOFC. 
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